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Sampling versus filtering in Large- Eddy 

simulations 

By 0. Debliquy, B. Knaepen & D. Carati f and A. A. Wrayj 


A LES formalism in which the filter operator is replaced by a sampling operator is 
proposed. The unknown quantities that appear in the LES equations originate only from 
inadequate resolution (discretization errors). The resulting viewpoint seems to make a 
link between finite difference approaches and finite element methods. Sampling opera- 
tors are shown to commute with nonlinearities and to be purely projective. Moreover, 
their use allows an unambiguous definition of the LES numerical grid. The price to pay 
is that sampling never commutes with spatial derivatives and the commutation errors 
must be modelled. It is shown that models for the discretization errors may be treated 
using the dynamic procedure. Preliminary results, using the Smagorinsky model, are very 
encouraging. 


1. Introduction 

The analysis presented here owes its origin to the unique practical motivation for 
Large-Eddy Simulations (LES’s), i.e. the lack of adequate computational power. Indeed, 
LES is a discretization of the Navier-Stokes equations on a grid which is too coarse. Tra- 
ditionally, the LES equations are presented as a smoothed version of the Navier-Stokes 
equations. The smoothing operation generates an unknown term that must be modelled. 
The solution of the LES equations also requires numerical schemes with, unavoidably, 
discretization errors. Since resolution is always the driving constraint in LES, the dis- 
cretization errors are expected to play a significative role (Ghosal 1996; Kravchenko & 
Moin 1997; Chow & Moin 2003). However, several attempts have been proposed in order 
to minimize their effect on the large scale dynamics (Vasilyev et al. 1998; Marsden et 
al. 2002). 

The viewpoint adopted in the approach developed here is very different, to such an 
extent that it might be viewed as provocative. This is however not the objective. The 
strategy that has been investigated is to abandon the idea of using smoothing filters. 
In that case, the discretization errors become the only term that requires a modelling 
effort. This viewpoint should not be confused with the MILES (Monotone Integrated 
Large Eddy Simulation) approach (Park et al. 2004) in which the numerical errors are 
expected to dissipate the correct amount of energy and to stabilize the under-resolved 
simulation. On the contrary, in the following, a real modelling effort is developed in 
order to account for the discretization errors. In particular, using the successful dynamic 
procedure, two embedded grids are used to optimize the model for the discretization 
errors. 
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In spectral codes, filtering or reducing the resolution are two indistinguishable oper- 
ations. In that sense, for these codes, the discretization errors and the modelling errors 
are identical and the methodology developed in this report would reduce to a simple 
re-interpretation of the LES equations without any practical consequence. For finite dif- 
ference codes the situation is different. The filter width A j and the grid spacing A g may 
be considered as two independent parameters. In the limit A g / Af — > 0, the discretiza- 
tion errors are supposed to be negligible with respect to the “sub-filter” term. Due to 
finite computational ressources, such a limit is however never used in practice. A rea- 
sonable compromise is to consider that A g / Af < 1. In the following, the opposite limit 
(A g / Af -* oo), obtained when no smoothing filter is used, is investigated for finite dif- 
ference schemes. A bridge between the this approach and the traditional viewpoint is also 
proposed. The discretization scheme is considered to result from a sampling operation 
that shares some properties with filters. 


2. Definitions 

The purpose of this section is to introduce the sampling operators and to make the 
link between two apparently incompatible viewpoints on LES. The first one is the tradi- 
tional filtering approach in which the Navier-Stokes velocity field u, is transformed into 
a continuous filtered velocity iq. The second viewpoint is the discrete representation u, 
of the velocity field on a numerical grid in the actual LES simulation. In the filtering 
approach, the LES velocity Ui is given by the following expression 

Ui(x) = J G(x,y) Ui(y) d 3 y, (2.1) 

in which, G(x, y) is the filter kernel. In general, the filtered field depends on values of the 
original velocity field from the entire integration domain (lei 3 . Very few constraints 
are imposed on the filter kernel G(x,y) except that G(x,y) = 0 for y # 0 and the 
normalization condition: 

J G(x,y)d 3 y = 1. (2.2) 

When considering the discrete representation of the signal u, on the LES grid, the sampled 
signal is not a function from ff into ffi anymore, but a mapping from the set of grid points 
{x e } into ffi. It is however possible to keep the same type of expression, 

Ui(x) = j H(x,y) Ui(y) d 3 y . (2.3) 

In order to make the link with sampling and finite difference schemes, the discrete-filtered 
function Ui(x) should depend only on the grid values of Ui(y), so that the general form 
of the kernel H has to be: 

H(x,y) = Y h l (x) S(y — x l ) , (2.4) 

i 

and the discrete-filtered signal is given by: 

Ui(x) = Y h l (x) Ui(x l ) . (2.5) 

i 

Such an expression amounts to assuming an interpolation of the signal between the grid 
points. However, only the values of the Navier-Stokes velocity on the grid are required 
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to construct U*. Imposing a normalization condition on the discrete-filtering operator is 
also desirable. Since only the grid point values of w, matter, it is sufficient to impose the 
normalization condition on the grid: 

f H(x l \y) d 3 y = 1, (2.6) 

J Li 

where the square matrix H ili2 = h l 2 (x tl ) has the dimensions N x N where N is the 
number of grid points. However, in general, the interpolation of a constant signal on a 
grid should also be a constant signal between the grid points and it is very natural to 
impose the normalization condition between the grid points as well: 

[ H(x,y) d 3 y =£V(f) = l. (2.7) 

J i 


The process of sampling a signal is very similar to a projection. However, the transfor- 
mation (2.5) does not necessarily correspond to the application of a projective operator 
to Ui(x). Indeed, if the signal is discrete-filtered again, the result is given by 

Uj{x) = T h ll {x) Yj h e2 (x (l ) Ui(x t2 ) , (2.8) 

h h 

and, introducing the IV-dimensional vectors = Ui(x il ), the discrete-filtering operator 
must satisfy the following constraint 

Y H tlt2 ul 2 = wf , (2.9) 

h 

in order to be considered as a projector. Remarkably enough, the same condition ensures 
that the discrete-filtered signal takes the same value as the Navier-Stokes signal on the 
grid points since the right-hand-side of relation (2.9) is simply ufix 12 ). The three vectors 
U } 2 ( i = 1,2,3) are thus eigenvectors of the square matrix H ilt2 with the eigenvalue 
A = 1. These constraints, for three independent signals u±, 112 and U 3 , are difficult to 
satisfy and, although no attempt will be done to prove it, the only solution seems to be: 

H llt 2 = 5 (li2 . (2.10) 


This choice automatically satisfies the normalization constraint as well. 

It is also interesting that, imposing the projector constraint, the derivatives of these 
functions at the grid points entirely define the numerical scheme that has to be used 
when solving the discrete-filtered Navier-Stokes equations. Indeed, the derivatives of Tii 
on the grid are given by 


d m Ui(x h ) = Dl m 2 ufix 12 ) = Y 
Li Li 

dmnMx* 1 ) = Y D ™n U ^) = Y 

Li Li 

where the differentiation matrices are given by 

dh l 


£)LiL 2 _ 


dx 
& 2 h l 


r\LiLi _ 
mn dx m dx, 


(■ 2 h )- 


D ^ 2 Ui(Z (2 ) , 

(2.11) 

Dl mn Ui(x l2 ) , 

(2.12) 


(2.13) 


(2.14) 



136 


0. Debliquy, B. Knaepen, D. Carati & A. A. Wray 
The normalization constraint imposes: 


I>™' 2=0 ’ 

(2.15) 

tl 


D ia 2 = o . 

/ y mn v • 

(2.16) 

£2 



These constraints are fairly easy to satisfy. 


3. Example 

The simplest example of a finite difference scheme is probably the second order centered 
difference scheme on a regular one dimensional grid. In that case, the grid points are given 
by x l = L l/N where l = 1 • • • N. The signals will be assumed to be periodic and two grid 
points are added to the grid: x° = x N and x N+1 = x 1 . The first and second derivatives 
of the discrete signal are given by: 

du(x e ) 
dx 

d 2 u(x t ) 
dx 2 

where A = L/N is the grid spacing. The choice of function h e (x) is not unique, even 
when imposing (3.1) and (3.2). In order to simplify the structure of the function h e (x), 
only compact support functions will be considered. Moreover, since a homogeneous grid 
is assumed, it is reasonable to assume that the functions h l have all the same shape and 
are characterized by a single function: 

h\x)=g(^) (3.3) 

Since the constraints (3.1) and (3.2) only concern the first two neighbor points, the 
support will be assumed to be [—3/2, 3/2] so that only three points are included in 
it. The connection points at y = ±3/2 might lead to discontinuities. However, as an 
additional constraint, the function g is assumed to be C 1 . The second derivative is thus 
allowed to be discontinuous, but only between grid points. In that case, it is actually 
much easier to define g(y) separately for three intervals as follows: 

9(y) = Pa(y ) ©(^ + y) ©(^ - y) 

+Pb(-y ) ©(^ + y) ®{-\-y)+Pb{y) ©(^-y) 0(“+y), (3.4) 

where Q(y) is the Heaviside function. This choice implies that g(y) = g(—y), which is 
reasonable for a centered difference scheme. The constraint (2.10) imposes that p o (0) = 1 
and pb( 1) = 0. Taking into account (1) the fact that the functions h ( are assumed to be 
all the same (3.3) and (2) the compact support of g, the normalization constraint reduces 
to g(y) + g(y + 1) + g(y — 1) = 1. Considering the structure (3.4), this imposes: 


u(x l+1 ) - u{xr *) _ aA _ 

2A = 6 x U ’ 

u(x t+1 ) - 2 u(x t ) +u(x t ~ 1 ) _ 
A2 = 


p a (y) +Pb(y + 1) +Pb(-y + 1) = l 


(3.5) 
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Figure 1. The function g(y) defined by the expressions (3. 4, 3. 6-3. 7). 


The lowest order polynomials that can be found to satisfy these constraints as well as 
the continuity of g and of its derivative at y = ±3/2 and y = ±1/2 are given by: 


p a {y) = 1 - y 2 - 4 y 4 , 

(3.6) 

Pb(y) =2(1+2/) (^ + y) 2 (2 + 5y+4y 2 ). 

(3.7) 


The resulting function g is shown in the Figure 1. For summarizing the properties of the 
discrete-filtering operator, the above choice for the function g ensures that (on the grid): 


u = u 
u 2 = u 2 , 

du(x e ) u(x l+1 ) -Ti^X 1 - 1 ) a A_, l\ 

—&T = 2A = dx U{X } ’ 

d 2 u(x ( ) _ u(x t+1 ) - 2u(x t ) +u(x ( - 1 ) _ A 2 ( 

Q x 2 ^2 W* ) u \ x ) 


and that (everywhere): 


u(x) = u(x ) . (3.8) 

For three-dimensional systems, the function h e (x) is simply a product of one-dimen- 
sional functions g and the grid point coordinate x e can be characterized by a vector of 
indices (l x ,l y ,l z )'- 


vM{x) = g ^ 
and x l = (x tx ,y iy ,z (z ). 


x — ar 


y-y 




(3.9) 


4. Interpolation and diagnostics 

In the preceding sections, interpolation functions compatible with a given finite differ- 
ence scheme have been introduced. They correspond to functions for which the deriva- 
tives on the grid are evaluated exactly by the spatial discretization scheme. The functions 
h/ 1 (x) thus play the same role for finite difference schemes as the plane waves e lk x for 
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Figure 2. The comparison between the sampling operator and a sharp Fourier cutoff is 
shown for a one dimensional signal. The solid line corresponds to the original periodic signal 
on the interval [0, 27r] built using 32 modes with random amplitudes. The dashed curve 
corresponds to the sampled signal on a 8 point grid reconstructed using the interpolation 
functions h l (x). The dotted line is obtained by applying a Fourier cutoff keeping 8 modes. 

spectral methods. In that sense, their introduction improves the mathematical ground- 
ing of finite differences, though they do not modify the practical implementation of 
these numerical methods. The approach developed here could then be considered at first 
sight as an interesting mathematical framework for finite difference schemes but without 
practical implications. There are however practical implications. The first one is in the 
post-processing of results obtained using a finite difference scheme. Many quantities are 
routinely computed in a numerical simulation. It is not the purpose of this work to derive 
the influence of the interpolation function on the exhaustive list of possible diagnostics. 
Only two important quantities, the energy and the dissipation will be considered. More- 
over, in order to simplify the discussion, the following developments are restricted to 
one-dimensional systems. 

By definition, the one-dimensional kinetic energy is given by 

E = ^ J dx u(x) u(x ) , (4.1) 

= Q ilk ( 4 - 2 ) 

Z ixll 

where 

Qiuh = j dx h ll (x ) h l2 (x). (4.3) 

In spectral methods, thanks to the Parseval theorem, the matrix Q tli2 reduces to 
It is possible to impose this constraint for finite difference schemes at the expense of an 
increase in the degree of the polynomials used in the definition of the functions h ll (i T). 
However, the choice for the set of most relevant diagnostics is problem dependent. Any 
attempt to simplify their expression would thus make the definition of the function h e (x) 
problem dependent as well. To avoid this undesirable feature, the “simplest” functions 
h (l (x) derived in the previous section and compatible with the continuity and normal- 
ization constraints as well as with the numerical scheme will be kept. Consequently, for 
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the one dimensional second order centered finite difference scheme, the matrix Q Cli2 is 
not diagonal but can be evaluated explicitly: 

Q* 1 ’* 2 =qi S h,k +q2 (^ 1,^+1 +q3 (6 (l ’ e 2+2 +<^-2) , (4.4) 

where qi = 11339/13860, qi = 121/1260 and q 3 = —47/9240. The property q\ + 2qi + 
2q 3 = 1 can be derived from the normalization properties of the functions h e (x). 
Similarly, the one dimensional dissipation is defined by: 

8 = 2v j dx { d x u(x )) ( d x u(x )) , (4.5) 

= 2 ^ u(x ei ) W llt2 u(x i2 ) , (4.6) 

hh 

where 

W 1 " 12 = J dx ( d x h (l (x )) ( d x h l2 (x )) . (4.7) 

The matrix W eit - 2 can also be evaluated and is given by: 

W* 1 -* 2 =w 1 6 tl ’ t2 +wi (6 (l ’ (2+1 + S* 1 ’* 2 - 1 ) + w 3 (5 ll ’ ta+2 +<5^’ <2 - 2 ) S: (4.8) 

where wi = 1009/315, W 2 = —179/105 and w 3 = 13/126. The property wi +2 W 2 +2 w 3 = 
0 can also be derived from the normalization properties of the functions h e (x). 


5. Dynamic Procedure 

Another place where the sampling viewpoint has a direct influence is in the implemen- 
tation of the dynamic procedure. Indeed, the definition of the velocity on the test-level 
grid requires the explicit definition of the coarsening operator. 

There is no reason to expect the dynamic Smagorinsky model to be a good candidate for 
representing the effect of the commutation error between the discrete-filtering operator 
and the spatial derivative. However, since this model has been thoroughly benchmarked 
and, at the very least, has the property to be dissipative so that it should stabilize the 
simulation, there is some rationale to evaluate and test the dynamic procedure for this 
model as a first step. 

Numerical simulations, both DNS and LES, can be seen as the time advancement of 
a discrete-filtered velocity. We will consider three grids: The DNS grid Q A °, the LES 
grid and a test grid Q^ 2 that is used only in the implementation of the dynamic 
procedure. In order to stress that the discrete-filtering operator is really based on the 
sampling concept, its application to the incompressible Navier-Stokes velocity is denoted 
by: 

u?°=S A °o Ui , (5.1) 

where uf° is the DNS field and the LES and test fields will be noted respectively uf 1 
and uf 2 . The three grids are assumed to be successively embedded Q A2 C f? Al C f/ A °. 
This is enough to guarantee that 

£A 2 = £A 2 0 £A! = £A 2 0 £A, q £A 0 ; (5.2) 

S Al = S Al O S A ° , (5.3) 

since only the grid point values of the velocity are used to construct the uf°, uf 1 and 
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up velocities. The DNS equation can be written as: 

d t up = -dpupup - SfV 0 + i/A A °u A ° i (5.4) 

-V ' 

Jvf°(« A o) 

where the differentiation operators <9 A ° and A A ° define the numerical scheme used in the 
DNS. By definition, the discretization error and, in particular, the difference between the 
discrete-filtered value of the velocity derivative (djUi) A ° and its evaluation on the grid 
dpup is supposed to be small enough so that it need not be modelled. This is obviously 
not the case in the LES equation: 

d t u p = -dpupup - dpp^ + i/A Al u Al +£ Al ’ A ° , (5.5) 

' v 

iV Al (« A i) 

and the additional term 


£p’*° = <S Al o lV A °(u A °) - Np (w Al ) , 


(5.6) 


has to be modelled. It contains three contributions from the pressure, the convective and 
the viscous terms. Actually, these terms should be constructed using the Navier-Stokes 
equations instead of the DNS equations. However, since it is assumed that the DNS 
velocity is a fair representation of the actual velocity, the expression (5.6) is acceptable 
and allows a more homogeneous presentation. Using the same notations, the test-level 
equations can be written as 



+£■ 




iVf 2 (« A 2) 


(5.7) 


and the additional term that has to be modelled at the test-level reads: 


£f 2 ’ A ° = S A2 o Np (u A ° ) - Np(u A2 ) . (5.8) 

Equivalently, the equations for up can be obtained by applying S* 2 on the LES equa- 
tions: 

d t up = S* 2 oNp(u^)+ 5 A2 o Sp ,A ° , (5.9) 

Consistency between these two formulations imposes the following identity: 


£»A 2 ,Ai 


£A 2 ,A() 


S^oSp'** , 


(5.10) 


where 

£p ’ Ai =<S A2 oNp(u^) -Np(u A2 ). (5.11) 

Equation (5.10) is nothing else the equivalent of the Germano identity for discrete- 
filtering operators, written at the vector level instead of the tensor level. The term (5.11) 
is known in terms of the LES velocity only (since up = S* 2 o up ) and plays the role of 
the Leonard tensor in the Germano identity. The right-hand-side depends on the terms 
that have to be modelled at both LES and test levels. 

If the Smagorinsky model is used at both LES and test levels, the identity is only an 
approximation that reads: 

£ A2 ’ Al w 2 C Mp (5.12) 


Mp = A 2 2 dp 


|S A2 | Sp -S^oAl dp |S* 


qAi 


where 


(5.13) 
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The tensor Sffi and Sffi represents the strain tensors evaluated on the LES and test 
grids respectively. The LES and test grid spacings are, respectively, Ai and A 2 . In this 
formulation, it is assumed that the Smagorinsky parameter C is constant in the entire 
volume and is the same at both LES and test levels. Its optimal value is given by: 


1 <£f 2 ’ Al Mf 2 ) 

2 <Mf 2 Mf 2 ) 


(5.14) 


Assuming that a = A 2 /Ai is a constant, the dynamic procedure can be used to compute 
directly the constant c = CAf : 


1 <ff 2 ' Al mf 2 ) 

2 {mf 2 mf 2 ) ’ 


where 


m- 


= a 2 df 2 |S a 


- S ^ o df 1 | S A 


13 


and the model that has to be implemented in the code reads: 


C A 1? A 0 


2 c <9 Al (I s* 




(5.15) 

(5.16) 

(5.17) 


6. Numerical results 


This model has been implemented in a pseudo-spectral code with modified wave vec- 
tors in order to mimic the behavior of a second-order centered finite difference scheme. 
Although this method is not computationally optimal, it has the advantage to be very 
easy to implement both for the DNS and for the LES. The computational domain is a 
cubic 2n x 2n x 2n box. The DNS grid, Q A ° , corresponds to a 256 3 resolution and the 

-»FD 

modified wave vectors for the finite difference scheme k are defined by: 


k FD = (fc FD , fc FD , fc FD ) = — (sin (kl A 0 ) , sin (k s y A 0 ) , sin (k s z A 0 )) , 


(. k z Y D = 


A 2 


sm 


kl A 0 


+ sin 


k s y A 0 


+ shr 


fc|A 0 


(6.1) 

(6.2) 


where k s is the usual wave vector in the pseudo-spectral code. Note that, in the finite 
difference scheme, the discretization of the second order derivative is not obtained as the 
“square” of the first order derivative and consequently different definitions are used for 
the modified wave vectors entering the first and the second order derivatives. The quantity 
(fe 2 ) FD is only used in the viscous term. The incompressibility condition is enforced by 
projecting the velocity w(k FD ) on the plane perpendicular to k FD . 

The sampling operators S Al and S A 2 are implemented straightforwardly in real space. 
It should be acknowledged however that enforcing incompressibility is an operation that 
now depends on the resolution. Consequently, it slightly modifies the properties of the 
sampling operators, and the velocity values on the grids G A °, Q Al and Q A 2 do not coincide 
exactly. The combination of the sampling and the incompressibility projections remains 
nevertheless a projective operator. 

The LES field is thus obtained by sampling using S Al on a 32 s grid. The computations 
in the implementation of the dynamic procedure that corresponds to the grids Q Al and 
Q A2 are then performed by a simple re-definition of the modified wave vectors. 

The preliminary tests, though made with a modified spectral code and a limited res- 
olution, are quite illuminating. First, it shows clearly that a model is required for the 
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Figure 3. Left : Energy decay. Comparison between the sampled DNS ( ), the LES 

curves without model ( ) and with the dynamic Smagorinsky model with the parameter 

C computed using both convective and viscous contributions ( ), the convective con- 
tribution only ( ) and the viscous contribution only ( ). Right : Energy spectrum 

at the end of the computation, same symbols, r is the large-eddy turnover time. 


decaying turbulence case. The under-resolved DNS (LES without model) exhibits an in- 
crease of the energy after a short decaying period. The dynamic Smagorinsky model has 
been implemented in three versions. Indeed, the “Leonard” vector (5.11) contains 

two contributions originating respectively from the convective-f-pressure terms and from 
the viscous term. In the computation of the dynamic Smagorinsky constant C, the purely 
convective, the purely viscous, and the full Leonard vectors have been implemented. Not 
surprisingly, the purely viscous version does not perform satisfactorily. It is actually even 
worse that the no-model case. Also, the purely convective version produces the best 
agreement with the sampled DNS. The final spectra are also presented. They all exhibit 
a piling-up of the energy in the high wave-number range. However, the no-model and the 
purely viscous dynamic model both overpredict the low wave-number energy range, while 
the dynamic models that account for the convective discretization errors both correctly 
predict the low k spectrum. 

Although it would be quite hazardous to make definitive and strong statements on the 
basis of these preliminary numerical results, they indicate that (1) the implementation 
of a dynamic procedure with sampling operators yields reasonable predictions in the 
homogeneous decaying turbulence test case and (2) the Smagorinsky model is acceptable 
for modelling the numerical discretization errors on the convective term but not for the 
viscous term. 


7. Discussion 

The methodology presented here is unusual in the sense that no filter is used to define 
the LES fields. The unknown quantities that appear in the LES equations thus originate 
only from the lack of resolution (discretization errors) and filters are replaced by sam- 
pling operators. The resulting viewpoint seems to make the link between finite difference 
approaches and finite element methods. The sampling operators have some very pleasant 
properties. They commute with nonlinearities, they are projections and they explicitly 
define the numerical grid. The obvious drawback is that they do not commute with spa- 
tial derivatives and the commutation errors must be modelled. It has been shown that 
models for the discretization errors may be treated using the dynamic procedure. Prelim- 
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inary results, using the Smagorinsky model, are very encouraging, though the modelling 
of discretization errors by an eddy viscosity term may be debatable. 

Considering the encouraging numerical results obtained in the very simple tests pre- 
sented in section VI, it can be concluded that the methodology in which a LES is defined 
by sampling operators deserves, at the very least, further investigations. It is however not 
the purpose of this last section to oversell the sampling-based LES approach. We have 
thus chosen to close this discussion by pointing out two difficulties that have been raised 
in the course of the Summer Program and that definitively deserve some attention. 

First, as mentioned before, the sampling operators do not preserve the incompressibil- 
ity of the flow. This difficulty has to be expected since the incompressibility condition 
involves spatial derivatives: 

d iUi = 0. (7.1) 

The non-commutation of the sampling operator and the spatial derivative thus gener- 
ates an unknown term in the sampled version of Eq. (7.1). This unknown term has been 
neglected in the preliminary tests presented here. The modelling of such a term would 
probably be fairly difficult and might change the nature of the numerical scheme it- 
self since the incompressibility condition is a very strong constraint in most numerical 
methods. To the best of our knowledge, the commutation error in the incompressibility 
condition has never been taken into account in more traditional LES in which non- 
commutative filters are also used. It might turn out that the only practical strategy is to 
ignore this effect, but numerical and theoretical investigations of the commutation error 
in the incompressibility condition would be welcome. 

The second difficulty has been pointed out by O. Vasilyev: the sampling operators do 
not conserve the mean of the original signal. This can easily be seen if a sine wave with 
zero mean is sampled systematically at the minima, the resulting sampled signal would 
be constant and negative. This is a special case of aliasing due to too coarse sampling. 
Whether this is an important issue or not is not yet clear. Obviously, it has no damaging 
consequences in the isotropic turbulence tests. However, the consequences in complex 
geometries with mean flows remain to be explored. 
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